Gamra: Simple meshing for complex earthquakes

The static offsets caused by earthquakes are well described by elastostatic models with a discontinuity in the displacement along the fault. A traditional approach to model this discontinuity is to align the numerical mesh with the fault and solve the equations using finite elements. However, this distorted mesh can be difficult to generate and update. We present a new numerical method, inspired by the Immersed Interface Method (Leveque and Li, 1994), for solving the elastostatic equations with embedded discontinuities. This method has been carefully designed so that it can be used on parallel machines on an adapted finite difference grid. We have implemented this method in Gamra, a new code for earth modeling. We demonstrate the correctness of the method with analytic tests, and we demonstrate its practical performance by solving a realistic earthquake model to extremely high precision. & 2016 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/). 1. Motivation A common feature of many earthquakes is a complex network of intersecting faults. Accurately modeling the static offsets and associated large scale deformation due to this fault geometry is crucial to a reliable understanding of seismic hazards (Marshall et al., 2008). The behavior of these faults is relatively well described by the equations of variable modulus elastostatics. However, for realistic faults, the displacement does not gradually taper off, but rather ends abruptly. This abrupt termination gives rise to a logarithmic singularity in the displacement (Okada, 1992). In realistic faults, these singularities are smoothed out by non-linear processes at the fault tips that are on a scale that are many orders of magnitude smaller than the fault itself. These characteristics make it challenging to numerically model realistic fault networks. In addition, elastostatics is only one piece of the puzzle when modeling the earthquake cycle. We want to incorporate an elastostatic solver into an overall algorithm for modeling the entire earthquake cycle (Barbot et al., 2012). We desire a unified method, using the same mesh, architecture, and boundaries, that can solve elliptic equations (for static offsets of earthquakes), parabolic equations (for poro-elastic and visco-elastic evolution between earthquakes), and hyperbolic equations (for dynamic rupture during an earthquake). Then we will have a powerful tool for self r Ltd. This is an open access articl y), consistent models of the entire earthquake cycle. At present, one relatively successful approach to building this kind of tool uses boundary integral methods (Barbot et al., 2012; Kaneko et al., 2010; Lapusta and Barbot, 2012; Hori et al., 2004; Kato, 2004; Matsuzawa et al., 2010; Rice, 1993; Shibazaki and Shimamoto, 2007; Smith and Sandwell, 2004). However, boundary integral methods inevitably make simplifications in the geometry or the physics of the problem. Finite-element methods (Aagaard et al., 2013; Hassani et al., 1997; Melosh and Williams, 1989; Puente et al., 2009; Kaneko et al., 2008, 2011) provide a natural way to fully represent the geometry and the physics as long as the mesh conforms to the faults. Generating these conforming meshes can be quite challenging and time consuming, especially when the faults intersect. The extended finite element method (Becker et al., 2009; Coon et al., 2011; Zangmeister, 2015) shows great promise in addressing this problem with mesh generation, though it has yet to be applied to realistic 3D earthquake models. Finite difference methods, on the other hand, have not traditionally been used for this kind of problem. Straightforward implementations of finite differences require that the displacement be continuous and differentiable. This limitation spurred the development of the Immersed Interface Method (IIM) (Leveque and Li, 1994). IIM explicitly models the discontinuous jump, resulting in a series of corrections to the ordinary finite difference stencils. IIM has spawned a number of variations, and some of these have been applied to various problems in elastostatics (Rutka et al., 2006; Rutka and Wiegmann, 2006; Botella and Cheny, 2010; Zhu et al., 2012). None of them have looked at models most relevant to e under the CC BY license (http://creativecommons.org/licenses/by/4.0/). W. Landry, S. Barbot / Computers & Geosciences 90 (2016) 49–63 50 earthquakes, where we prescribe the discontinuity in the displacement. More importantly, none of them have discussed how to handle the difficulties associated with the singularity at the fault tip. Finally, none of these methods have been implemented on adapted grids or parallel machines. The purpose of this paper is to describe a new method, inspired by IIM, that naturally handles all of the difficulties associated with faults. This method was developed with an eye towards performance, so it naturally extends to the use of parallel machines and highly adapted grids. With this solver in place, we can then use the existing deep understanding of how to implement hyperbolic and parabolic solvers for the equations specific to earthquakes in a finite difference framework (Day, 1982; Dunham and Archuleta, 2005; Dunham et al., 2011; Andrews, 2002; Day et al., 2005; Harris et al., 2009; Olsen et al., 1997; Ely et al., 2009, 2010; Cui et al., 2010; Kozdon et al., 2013; Moczo et al., 2014). We first describe the equations of linear elasticity, how we treat internal dislocations, and how we solve these equations on an adapted mesh. Then we demonstrate the correctness of the method and our implementation with a series of analytic tests. Finally, we document the performance of our implementation with a simulation of the 1992 Mw 7.3 Landers earthquake. The algorithm described in this paper is implemented in Gamra, a code available at https://bitbucket.org/wlandry/gamra. Gamra is a French acronym for Géodynamique Avec Maille Rafinée Adaptivement, meaning “geodynamics with adaptive mesh refinement”. Fig. 1. Reference cell showing where the displacement and moduli are defined. The bottom left is at x1⁄40, y1⁄40, and the top right is at δ = x x , δ = y y. 2. Methods We begin by describing the equations of linear elasticity (Section 2.1) and the mesh we use for solving them (Section 2.2). Then we describe the Gauss–Seidel smoother that we use as a component in our solvers (Section 2.3). Then we describe the corrections we make to treat internal dislocations of arbitrary orientation in two and three dimensions (Section 2.4). Then we describe how we implement boundary conditions (Section 2.5). With these components, we have a stable, accurate solver for earthquake physics. However, this will not be a fast solver without multigrid. To implement multigrid (Section 2.6), we need coarsening (Section 2.6.1) and refinement (Section 2.6.2) operators. To implement adaptive multigrid, we also need to set boundary conditions at coarse-fine boundaries (Section 2.6.3). 2.1. Governing equations We solve the Navier's equation for elastostatic deformation with the infinitesimal strain approximation σ + = ( ) f 0, 1 ji j i , where the stress components sji are defined using Hooke's law in terms of the displacement components vi, Lame's first parameter λ, and the shear modulus μ σ μ δ λ (→) ≡ ( + ) + ( ) v v v v . 2 ji i j j i ij k k , , , We use Einstein summation notation, where each index i, j, k is understood to x, y, and z in turn, repeated indices are summed, and commas (,) denote derivatives. For all of our test problems, the stress tensor will be symmetric σ σ ( = ) ij ji . In addition, the forcing term fi is zero for many of our test problems. But equivalent body forces can be used to represent inelastic deformation in quasi-static deformation simulations (Barbot and Fialko, 2010; Rousset et al., 2015; Rollins et al., 2015). Therefore the inclusion of body forces in Eq. (1) is critical for modeling quasi-static deformation due to off-fault processes. 2.2. Staggered grid We discretize the equations on a staggered grid, with the displacement located at cell faces as shown in Fig. 1. Our method requires the shear modulus (μ) at both the cell centers and cell corners. Since μ is a given function of space, we could compute it exactly at both cell centers and corners. We have found that we get larger reductions in the residuals for each multigrid V-cycle by using the given function to compute the cell centers, and then using the geometric mean to fill the value at the cell corners. Specifically, in 2D, for a reference cell where the bottom left corner is located at x1⁄40, y1⁄40, μ at that corner is ( ) μ μ μ μ μ = ( ) δ δ δ δ δ δ δ δ − − − − . 3 x y x y x y x y 0,0 /2, /2 /2, /2 /2, /2 /2, /2 1/4 The subscripts δy 0, /2 indicate the variable located at an offset of x1⁄40, δ = y y/2 from the bottom left corner. So 0,0 is the bottom left corner, δy 0, /2 is the left face, and δ δ x y /2, /2 is the cell center. The Lame parameter λ is only needed at cell centers, so there is no extra interpolation step. We can specify μ and λ one of two ways: analytic expressions and tables. We use the muparser library (Berg, 2014) to evaluate analytic expressions. To compute the modulus at the boundary, we may need the modulus at a point outside the boundary. For analytic expressions, we evaluate the expression at that outside point. For moduli given by a table, we choose the closest point covered by the table. For multigrid, the modulus on coarser levels is interpolated from finer levels, not directly computed. Using the interpolated values rather than the directly computed values results in larger reductions in the residuals for each multigrid V-cycle. The interpolation onto the cell centered modulus is a simple arithmetic average of all of the fine points in the coarse cell. This treatment of the modulus works well for the moderate jumps in material properties seen in realistic models of Fig. 2. Update schedule for Gauss–Seidel relaxation in 2D. Updates for 3D follow a similar pattern. W. Landry, S. Barbot / Computers & Geosciences 90 (2016) 49–63 51 earthquake regions. More extreme jumps would require a more sophisticated treatment, such as applying IIM to material interfaces as well as faults. 2.3. Gauss–Seidel relaxation The core of the solver is a red-black Gauss-Seidel relaxation. We first define the residual as the non-zero remnant of Eq. (1): σ (→ → ) = + ( ) r v f f , . 4 i ji j i , We discretize the residual in the usual way with centered differences. To be explicit, in 2D, we write the x component as ( ) ( ) ( ) ( ) σ λ μ λ μ = + + + ( + ) v v v v 2 . jx j x x x y y x x y y x y , , , , , , , , where, in the reference cell ⎡ ⎣⎢ ⎤ ⎦⎥ ( )( ) ( )( ) ( ) ( ) λ μ